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Abstract Over the past years, various representation systems which sparsely ap- 
proximate functions governed by anisotropic features such as edges in images have 
been proposed. We exemplarily mention the systems of contourlets, curvelets, and 
shearlets. Alongside the theoretical development of these systems, algorithmic re- 
alizations of the associated transforms were provided. However, one of the most 
common shortcomings of these frameworks is the lack of providing a unified treat- 
ment of the continuum and digital world, i.e., allowing a digital theory to be a natural 
digitization of the continuum theory. In fact, shearlet systems are the only systems 
so far which satisfy this property, yet still deliver optimally sparse approximations 
of cartoon-like images. In this chapter, we provide an introduction to digital shearlet 
theory with a particular focus on a unified treatment of the continuum and digital 
realm. In our survey we will present the implementations of two shearlet transforms, 
one based on band-limited shearlets and the other based on compactly supported 
shearlets. We will moreover discuss various quantitative measures, which allow an 
objective comparison with other directional transforms and an objective tuning of 
parameters. The codes for both presented transforms as well as the framework for 
quantifying performance are provided in the Matlab toolbox [ShearLabj 



1 Introduction 

One key property of wavelets, which enabled their success as a universal method- 
ology for signal processing, is the unified treatment of the continuum and digital 
world. In fact, the wavelet transform can be implemented by a natural digitization of 
the continuum theory, thus providing a theoretical foundation for the digital trans- 
form. Lately, it was observed that wavelets are however suboptimal when sparse 
approximations of 2D functions are seeked. The reason is that these functions are 
typically governed by anisotropic features such as edges in images or evolving shock 
fronts in solutions of transport equations. However, Besov models - which wavelets 
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optimally encode - are clearly deficient to capture these features. Within the model 
of cartoon-like images, introduced by Donoho in |fT3l in 1999, the suboptimal be- 
havior of wavelets for such 2D functions was made mathematically precise; see also 
Chapter |4|. 

Among the various directional representation systems which have since then 
been proposed such as contourlets [12|, curvelets |9|, and shearlets, the shearlet 
system is in fact the only one which delivers optimally sparse approximations of 
cartoon-like images and still also allows for a unified treatment of the continuum 
and digital world. One main reason in comparison to the other two mentioned sys- 
tems is the fact that shearlets are affine systems, thereby enabling an extensive theo- 
retical framework, but parameterize directions by slope (in contrast to angles) which 
greatly supports treating the digital setting. As a thought experiment just note that 
a shear matrix leaves the digital grid invariant, which is in general not true for 
rotation. 

This raises the following questions, which we will answer in this chapter: 

(PI) What are the main desiderata for a digital shearlet theory? 
(P2) WTiich approaches do exist to derive a natural digitization of the continuum 
shearlet theory? 

(P3) How can we measure the accuracy to which the desiderata from (PI) are 
matched? 

(P4) Can we even introduce a framework within which different directional trans- 
forms can be objectively compared? 

Before delving into a detailed discussion, let us first contemplate about these 
questions on a more intuitive level. 



1.1 A Unified Framework for the Continuum and Digital World 

Several desiderata come to one's mind, which guarantee a unified framework for 
both the continuum and digital world, and provide an answer to (PI). The following 
are the choices of desiderata which were considered in [20., 14.1 : 

• Parseval Frame Property. The transform shall ideally have the tight frame prop- 
erty, which enables taking the adjoint as inverse transform. This property can be 
broken into the following two parts, which most, but not all, transforms admit: 
o Algebraic Exactness. The transform should be based on a theory for digital 

data in the sense that the analyzing functions should be an exact digitization 
of the continuum domain analyzing elements, 
o Isometry of Pseudo-Polar Fourier Transform. If the image is first mapped 
into a different domain - here the pseudo-polar domain -, then this map 
should be an isometry. 

• Space-Frequency-Localization. The analyzing elements of the associated trans- 
form should ideally be highly localized in space and frequency - to the extent to 
which uncertainty principles allow this. 
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• True Shear Invariance. Shearing naturally occurs in digital imaging, and it can 
- in contrast to rotation - be precisely realized in the digital domain. Thus the 
transform should be shear invariant, i.e., a shearing of the input image should be 
mirrored in a simple shift of the transform coefficients. 

• Speed. The transform should admit an algorithm of order 0{N^\ogN) flops, 
where A^^ is the number of digital points of the input image. 

• Geometric Exactness. The transform should preserve geometric properties par- 
allel to those of the continuum theory, for example, edges should be mapped to 
edges in transform domain. 

• Robustness. The transform should be resilient against impacts such as (hard) 
thresholding. 



1.2 Band-Limited Versus Compactly Supported Shearlet 
Transforms 

In general, two different types of shearlet systems are utilized today: Band-limited 
shearlet systems and compactly supported shearlet systems (see also Chapters [1] 
and ID). Regarding those from an algorithmic viewpoint, both have their particular 
advantages and disadvantages: 

Algorithmic realizations of the band-limited shearlet transform have on the one 
hand typically a higher computational complexity due to the fact that the window- 
ing takes place in frequency domain. However, on the other hand, they do allow a 
high localization in frequency domain which is important, for instance, for handling 
seismic data. Even more, band-limited shearlets do admit a precise digitization of 
the continuum theory. 

In contrast to this, algorithmic realizations of the compactly supported shearlet 
transform are much faster and have the advantage of achieving a high accuracy in 
spatial domain. But for a precise digitization one has to lower one's sights slightly. 
A more comprehensive answer to (P2) will be provided in the sequel of this chapter, 
where we will present the digital transform based on band-limited shearlets intro- 
duced in fSOl and the digital transform based on compactly supported shearlets from 
I22J. 



1.3 Related Work 

Since the introduction of directional representation systems by many pioneer re- 
searchers (||8]|9l[l0l|TTl[T2l), various numerical implementations of their directional 
representation systems have been proposed. Let us next briefly survey the main fea- 
tures of the two closest to shearlets, which are the contourlet and curvelet algo- 
rithms. 
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• Curvelets [7|. The discrete curve let transform is implemented in the software 
package CurveLab, which comprises two different approaches. One is based on 
unequispaced FFTs, which are used to interpolate the function in the frequency 
domain on different tiles with respect to different orientations of curvelets. The 
other is based on frequency wrapping, which wraps each subband indexed by 
scale and angle into a fixed rectangle around the origin. Both approaches can 
be realized efficiently in 0{N^logN) flops with being the image size. The 
disadvantage of this approach is the lack of an associated continuum domain 
theory. 

• Contourlets lfT2l . The implementation of contourlets is based on a directional fil- 
ter bank, which produces a directional frequency partitioning similar to the one 
generated by curvelets. The main advantage of this approach is that it allows a 
tree-structured filter bank implementation, in which aliasing due to subsampling 
is allowed to exist. Consequently, one can achieve great efficiency in terms of 
redundancy and good spatial localization. A drawback of this approach is that 
various artifacts are introduced and that an associated continuum domain theory 
is missing. 

Summarizing, all the above implementations of directional representation systems 
have their own advantages and disadvantages; one of the most common shortcom- 
ings is the lack of providing a unified treatment of the continuum and digital world. 

Besides the shearlet implementations we will present in this chapter, we would 
like to refer to Chapter |2| for a discussion of the algorithm in |16| based on the 
Laplacian pyramid scheme and directional filtering. It should be though noted that 
this implementation is not focussed on a natural digitization of the continuum theory 
and that the code was not made publicly available, both of which are crucial aspects 
of the work presented in the sequel. We further would like to draw the reader's 
attention to Chapter |3| which is based on ||2TI aiming at introducing a shearlet 
MRA from a subdivision perspective. Finally, we should mention that a different 
approach to a shearlet MRA was recently undertaken in ifTTl . 



1.4 Framework for Quantifying Performance 

A major problem with many computation-based results in applied mathematics is 
the non-availability of an accompanying code, and the lack of a fair and objective 
comparison with other approaches. The first problem can be overcome by following 
the philosophy of 'reproducible research' ifTSll and making the code publicly avail- 
able with sufficient documentation. In this spirit, the shearlet transforms presented 
in this chapter are all downloadable from http : //www . shearlab . org One 
approach to overcome the second obstacle is the provision of a carefully selected 
set of prescribed performance measures aiming to prohibit a biased comparison on 
isolated tasks such as denoising and compression of specific standard images like 
'Lena', 'Barbara', etc. It seems far better from an intellectual viewpoint to care- 
fully decompose performance according to a more insightful array of tests, each 
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one motivated by a particular well-understood property we are trying to obtain. In 
this chapter we will present such a framework for quantifying performance specifi- 
cally of implementations of directional transforms, which was originally introduced 
in II20I [T4II . We would like to emphasize that such a framework does not only pro- 
vide the possibility of a fair and thorough comparison, but also enables the tuning 
of the parameters of an algorithm in a rational way, thereby providing an answer to 
both (P3) and (P4). 



1.5 ShearLab 

Following the philosophy of the previously detailed thoughts, IshearLatjfJ was in- 
troduced by Donoho, Shahram, and the authors. This software package contains 

• An algorithm based on band-limited shearlets introduced in |20|. 

• An algorithm based on compactly supported shearlets introduced in 1221 . 

• A comprehensive framework for quantifying performance of directional repre- 
sentations in general. 

This chapter is also devoted to provide an introduction to and discuss the mathemat- 
ical foundation of these components. 



1.6 Outline 

In Section [2j we introduce and analyze the fast digital shearlet transform FDST, 
which is based on band-limited shearlets. Section[3]is then devoted to the presenta- 
tion and discussion of the digital separable shearlet transform DSST and the digital 
non-separable shearlet transform DNST. The framework of performance measures 
for parabolic scaling based transforms is provided in Section]?] In the same sec- 
tion, we further discuss these measures for the special cases of the three previously 
introduced transforms. 



2 Digital Shearlet Transform using Band-Limited Shearlets 

The first algorithmic realization of a digital shearlet transform we will present, 
coined Fast Digital Shearlet Transform (FDST), is based on band-limited shear- 
lets. Let us start by defining the class of shearlet systems we are interested in. Re- 
ferring to Chapter 11], we will consider the cone-adapted discrete shearlet system 
SH{(j),\if,\lr;A,A,A)^<P{(j);A)Ul'{\lf;A)U'P{f;A) with 4 = and 
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A=A = {U,k,m) : ; >0,|;t| < 2',ot(E Z^}. 

We wish to emphasize that this choice relates to a scaUng by 4^ yielding an integer 
valued parabolic scaling matrix, which is better adapted to the digital setting than a 
scaling by We further let be a classical shearlet (iff likewise with ,^2) — 
Wi^2,^i)), i.e., 

VA(^) = V^(^i,fe) = V/i(^i)V^2(|), (1) 

where Xj/i e L^{R) is a wavelet with e C°°(M) and supp C [-4,-i] U [5, 4], 
and V^2 S L?{M.) a 'bump' function satisfying 1/^2 € C"(M) and supp v>2 C [-1, 1]. We 
remark that the chosen support deviates slightly from the choice in the introduction, 
which is however just a minor adaption again to prepare for the digitization. Further, 
recall the definition of the cones "^n - ^^22 from Chapter [ 1 1. 

The digitization of the associated discrete shearlet transform will be performed 
in the frequency domain. Focussing, on the cone '^#21, say, the discrete shearlet trans- 
form is of the form 

/ ^ (/, V/,) = (/, = (/,2-^5 xif{SlA,-j-)eM^^-^'^"-)), (2) 

where rj = (j,k,m,i) indexes scale j, orientation k, position m, and cone i. Con- 
sidering this shearlet transform for continuum domain data (taking all cones into 
account) implicitly induces a trapezoidal tiling of frequency space which is evi- 
dently not cartesian. A digital grid perfectly adapted to this situation is the so-called 
'pseudo-polar grid', which we will introduce and discuss subsequently in detail. Let 
us for now mention that this viewpoint enables representation of the discrete shearlet 
transform as a cascade of three steps: 

1) Classical Fourier transformation and change of variables to pseudo-polar coor- 
dinates. 

2) Weighting by a radial 'density compensation' factor 

3) Decomposition into rectangular tiles and inverse Fourier transform of each tiles. 
Before discussing these steps in detail, let us give an overview of how these steps 



will be faithfully digitized. First, it will be shown in Subsection 2.1 that the two 
operations in Step 1) can be combined to the so-called pseudo-polar Fourier trans- 
form. An oversampling in radial direction of the pseudo-polar grid, on which the 
pseudo-polar Fourier transform is computed, will then enable the design of 'density- 
compensation-style' weights on those grid points leading to Steps 1) & 2) being an 



isometry. This will be discussed in Subsection 2.2 Subsection 2.3 is then concerned 
with the digitization of the discrete shearlets to subband windows. Notice that a 
digital analog of (|2]) moreover requires an additional 2D-iFFT. Thus, concluding 
the digitization of the discrete shearlet transform will cascade the following steps, 
which is the exact analogy of the continuum domain shearlet transform (|2]l: 

(SI) PPFT: Pseudo-polar Fourier transform with oversampling factor in the radial 
direction. 
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(52) Weighting: Multiplication by 'density-compensation-style' weights. 

(53) Windowing: Decomposing the pseudo-polar grid into rectangular subband 
windows with additional 2D-iFFT. 

With a careful choice of the weights and subband windows, this transform is an 
isometry. Then the inverse transform can be computed by merely taking the adjoint 
in each step. A final discussion on the FDST will be presented in Subsection|2.4| 



2.1 Pseudo-Polar Fourier Transform 

We start by discussing Step (SI). 

2.1.1 Pseudo-Polar Grids with Oversampling 

In 01, a fast pseudo-polar Fourier transform (PPFT) which evaluates the discrete 
Fourier transform at points on a trapezoidal grid in frequency space, the so-called 
pseudo-polar grid, was already developed. However, the direct use of the PPFT is 
problematic, since it is - as defined in 151 - not an isometry. The main obstacle 
is the highly nonuniform arrangement of the points on the pseudo-polar grid. This 
intuitively suggests to downweight points in regions of very high density by using 
weights which correspond roughly to the density compensation weights underly- 
ing the continuous change of variables. This will be enabled by a sufficient radial 
oversampling of the pseudo-polar grid. 

This new pseudo-polar grid, which we will denote in the sequel by Qr to indicate 
the oversampling rate R, is defined by 

QR = Q\unl, (3) 

where 

rt\ rr In U ln\ . N ^ fj ^ N RN ^ ^ RN\ /^x 

r)2 r/2n In 2i\ . N ^ n ^ N RN ^ „ ^ RN \ 

= 7v) : "2 2' ^" ^ T"}- (5) 

This grid is illustrated in Fig. [T] We remark that the pseudo-polar grid introduced 
in |5| coincides with Qr for the particular choice /? = 2. It should be emphasized 
that Qr ^ QgU Qg is not a disjoint partitioning, nor is the mapping (n, i?) i-> (— ^ • 
^7x)or(^,-^-^) injective. In fact, the center 

■^ = {(0,0)} (6) 

appears + 1 times in as well as i2|, and the points on the seam lines 
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Fig. 1 The pseudo-polai- grid Qr = Qj^U 0,1 for Af = 4 and = 4. 

=^i = {(|,-|):-f <«<f,«^0}. 
appear in both ^2^ and Q,^. 

Definition 1. Let N,R be positive integer, and let be the pseudo-polar grid given 
by ([3]). For an N x N image / := {/(m, v) :— y<M,v<^ — 1}, the pseudo-polar 
Fourier transform (PFFT) I of I evaluated on Qr is then defined to be 

ii,v=-N/2 

where mo> N is an integer. 

We wish to mention that mq >N is typically set to be mq — | {RN + 1 ) for com- 
putational reasons (see also Q), but we for now allow more generality. 

2.1.2 FastPPFT 

It was shown in |5 1, that the PPFT can be reaHzed in 0{N^ log A^) flops with N xN 
being the size of the input image. We will now discuss how the extended pseudo- 
polar Fourier transform as defined in Definition [T] can be computed with similar 
complexity. 

For this, let / be an image of size N x N. Also, mo is set - but not restricted - 
to be mo = ^{RN + 1); we will elaborate on this choice at the end of this subsec- 
tion. We now focus on £2^, and mention that the PPFT on the other cone can be 
computed similarly. Rewriting the pseudo-polar Fourier transform from Definition 
[l] for («i,«2) = (-1 • I, I) e ^2^, we obtain 
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-2n 



u.v=-N/2 

N/2-1 N/2-1 2.,, 
M=-Af/2v=-A'/2 

7V/2-1 / N/2-1 \ 

= L L /(M,v)e-lTl (7) 

u=-N/2 \v=-N/2 ) 

This rewritten form, i.e., (|7]l, suggests that the pseudo-polar Fourier transform / 
of / on Q.\ can be obtained by performing the ID FFT on the extension of / along 
direction v and then applying a fractional Fourier transform (frFT) along direction 
u. To be more specific, we require the following operations: 

Fractional Fourier Transform. For c £ C'^+', the (unaliased) discrete fractional 
Fourier transform by a £ Cis defined to be 



N/2 



J=-N/2 

It was shown in ||6l, that the fractional Fourier transform F^^^^c can be computed 
using 0{NlogN) operations. For the special case of a = 1 /{N + 1), the fractional 
Fourier transform becomes the (unaliased) ID discrete Fourier Transform (ID FFT), 
which in the sequel will be denoted by Fi . Similarly, the 2D discrete Fourier Trans- 
form (2D FFT) will be denoted by F2, and the inverse of the F2 by ' (2D iFFT). 

Padding Operator. For even, m > N an odd integer, and c G C'^, the padding 
operator E„^„ gives a symmetrically zero padding version of c in the sense that 



c[k) k — 2 J ■ ■ • I 2 ' 





2 ' ■ ■ ■ ' 2 

Using these operators, (|7]l can be computed by 

N/2-1 

/((»i,a)2) = £ FioERfi+i,NoI{u,n)e 

u=-N/2 



— 271U 



(RN+l)-N/2 



N/2 _2„ 



~N/2 

CCfi \ 

'N+V 



= {F""Ji;n)m, 



where / = £w+i,woFi o£k^+i,^o/ e C(«^+')><(^+i) and a„ - - (^^+"1)^/2 ■ Since 
the ID FFT and ID frFT require only 0{NlogN) operations for a vector of size 
A^, the total complexity of this algorithm for computing the pseudo-polar Fourier 
transform from Definition [I] is indeed 0{N^\ogN) for an image of size N xN. 
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We would like to also remark that for a different choice of constant mo, one can 
compute the pseudo-polar Fourier transform also with complexity 0{N^\ogN) for 
an image of size N xN. This however requires application of the fractional Fourier 
transform in both directions u and v of the image, which results in a larger constant 
for the computational cost; see also [6 |. 



2.2 Density-Compensation Weights 

Next we tackle Step (S2), which is more delicate than it might seem, since the 
weights will not be derivable from simple density compensation arguments. 

2.2.1 A Plancherel Theorem for the PPFT 

For this, we now aim to choose weights w : £2r IR+ so that the extended PPFT 
from Definition[T]becomes an isometry, i.e., 

^72-1 

£ |/(m,v)P= X W(«i,0>2)-|/>1,(02)|2. (8) 

u,v=-N/2 (mi ,(02)eiiR 

Observing the symmetry of the pseudo-polar grid, it seems natural to select weight 
functions w which have full axis symmetry properties, i.e., for all (0)1,0)2) G ^2^, 
we require 

w(0)i,0)2) =w(0>2,0)i), w(0)i,0)2) = w(-0Ji,0J2), w(0)i , 0)2) = w(0)i , -0)2). (9) 

Then the following 'Plancherel theorem' for the pseudo-polar Fourier transform on 
Qr - similar to the one for the Fourier transform on the cartesian grid - can be 
proved. 

Theorem 1 Let N be even, and let w: Qr ^ be a weight function satis- 

fying Q. Then ([8]) holds, if and only if, the weight function w satisfies 

5{u,v) = w(0,0) 

RN/2 

+ 4- L i: w(f,|.-^).cos(27r«.4;^).cos(27rv.4^.f) 

(=0.N/2 n=\ 
N/2-IRN/2 

+ 8- E j:wef,?i-^)-cosi27:u-^)-cosi27tv^-f) m 
forall -N+l < u,v<N-L 

Proof. We start by computing the right hand side of ([8|: 
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^ w((0i,(02)-|/(«l, 0)2)1^ 



^ w(a)i,a)2)- 

((Bl,(B2)Gi2R 

w(a)i,o)2)- 

((Bl,(B2)Gi2R 



£ /(«,V)/(«',V')' 



iV/2-l , . 

I{u,v)e "'0^ 

u.y=-N/2 

N/l-l N/2-1 

^ ^ I{u,v)I{u' y)e "0 

.H,v=-Ar/2H',v'=-A'/2 
^72-1 

= £ w(«i,0)2)- L 

. . -^(("-« )0'l + (v-v')<B2) 



u.r.ii'.\'' = -N/2 
{u,v)^{u' ,V') 



Choosing / = c„|,,,j 5(m — mi, v — vi) +c„2.v2^(m — M2, v — V2) for all —N/2 < mi , v'l, 
M2, V2 < N/2 — 1 and for all c„j,vj ,c„2,v2 G C!, we can conclude that ^ holds if and 
only if 



(fl)i,C02)Gi3;; 



5(m,v), <m,v<A^-1. 



By the symmetry of the weights (|9]), this is equivalent to 

Y w(a)i,a)2)-[cos(^MOJi)cos(^va)2)] = 5(m,v) 
(ffli,(B2)ei2R 



(11) 



for all —N+ 1 <M,v<A^— 1. From this, we can deduce that ( fTTj l is equivalent to 
( [Tol l, which proves the theorem. □ 

Notice that ([TOjl is a linear system with RN^/4 + RN/2 + 1 unknowns and (2N — 
1)^ equations, wherefore, in general, one needs the oversampling factor R to be at 
least 16 to enforce solvability. 



2.2.2 Relaxed Form of Weight Functions 



The computation of the weights satisfying Theorem [T] by solving the full linear 
system of equations ([T0| is much too complex. Hence, we relax the requirement 
for exact isometric weighting, and represent the weights in terms of undercomplete 
basis functions on the pseudo-polar grid. 

More precisely, we first choose a set of basis functions wi, . . . ,w„g : Qr — > 
such that 

"0 

Ywj{couCO2)^0 forall (tOi,t(}2) G X2«. 
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We then represent weight functions w : i2« — ^ M+ by 

«o 

w:=^C;Wy, (12) 

7=1 

with ci,...,c„(j being nonnegative constants. This approach now enables solving 
( [Tol l for the constants ci , . . . , c„q using the least squares method, thereby reducing the 
computational complexity significantly. The 'full' weight function w is then given 
by ([T2|. 

We next present two different choices of weights which were derived by this 
relaxed approach. Notice that {coi , (Oi) and {n,£) will be used interchangeably. 

Choice 1. The set of basis functions wi , . . . , is defined as follows: 
Center. 

= 1(0,0) I 

Boundary: 

^2 = l{((Ui,(U2):|«|=A'fi/2,(Ui=(B2} W3 = 1 {(a)i ,<»2):|«|=A'«/2, £1)17^61)2} i 

Seam lines: 

W4 — \n\ ■ l{((B,^(D2):i<|„|<A?«/2,(Oi=(i>2}i 

Interior: 

W5 = |n| • l{((0,,(02):l<l«l<A'fi/2,Oi/(i>2}- 

Choice 2. The set of basis functions wi , . . . ,Wi^ ^2+2 is defined as follows: 
Center: 

W'l = 1(0,0) I 

Radial Lines: 

Wt+2 ^ l{((Ui^(i>,):l<|„|<iVff/2,fi,2=^(Bi}' ^ = 0,l,...,A^/2. 

The associated weight functions are displayed in Fig. |2] In general, suitable 
weight functions usually obey the pattern of linearly increasing values along the 
radial direction. Thus, this is a natural requirement for the basis functions. 



2.2.3 Comparison of Weights 

A visual comparison shows that the patterns of the weight functions associated with 
Choices 1 and 2 are seemingly similar (see Fig.|2]i. However, carefully chosen mea- 
sures reveal that their performances can in fact be quite different. 

One essential criterion for the quality of a weight function is the degree to which 
it allows a Plancherel theorem for the pseudo-polar Fourier transform as studied in 
Theorem [T| This can be measured in the following way - the reader might want to 



compare this performance measure with the measures introduced in Subsection 4.2 
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Choice 1 Choice 2 



Fig. 2 Weight functions on the pseudo-polar grid for N = 256 and S = 8. 

Let P and P* denote the operators for the pseudo-polar Fourier transform and its 
adjoint, respectively, and let w - by slightly abusing notation - denote the weighting 
operator on the pseudo-polar grid ^2^;. Letting /? = 8, a sequence of 5 random images 
/i, . . ., /j of size N X N with standard normally distributed entries is generated to 
compute 

The performance of the weight functions arising from Choices 1 and 2 with respect 
to this measure is presented in Table [T] 



Table 1 Comparison of Choices 1 and 2 based on performance measure Mi . 



N 


32 


64 


128 


256 


512 


Choice 1 


4.2E-3 


4.0E-3 


1.8E-3 


1.5E-3 


8.8E-4 


Choice 2 


9.8E-3 


6.2E-3 


3.4E-3 


2.1E-3 


N/A 



Interestingly, a structured image, e.g., by using the measure 

||P*wP/-/||2 . . , , 

M2 := jT-j] , /the image Barbara , 

l|/||2 

yields an even better performance and a better distinction, see Table |2] One could 
reason that this behavior is due to the fact that the energy of most 'real' images 
is concentrated in the low frequency region, in which density compensation of the 
pseudo-polar grid is not as necessary as in the high frequency regions. 

These two tables show firstly, that with growing A^, the weighted pseudo-polar 
Fourier transform seems to converge to being an isometry on the testing image 
class. Secondly, judging from the relatively small deviation from being an isometry, 
it seems quite reasonable to choose a basis of weight functions forcing the weights 
to linearly increase along the radial direction. And, thirdly, although Choice 2 con- 
tains many more basis functions than Choice 1, the performance results are worse. 
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Table 2 Comparison of Choices 1 and 2 based on the performance measure M2 . 



N 


32 


64 


128 


256 


512 


Choice 1 


2.8E-3 


1.2E-3 


8.3E-4 


3.9E-4 


1.5E-4 


Choice 2 


5.6E-3 


2.8E-3 


2.2E-3 


9.1E-4 


N/A 



which is very counterintuitive. The reason for this is the numerical instability when 
computing a minimizing set of coefficients for the basis of weight functions, which 
causes these effects. 



2.2.4 Computation of the Weighting 



For the FDST - as also in the implementation in ShearLab - the coefficients 
in the expansion ( [T2| ) will be computed off-line, and then hardwired in the code. 
This enables the weighting of a function on the pseudo-polar grid to simply be a 
point- wise multiplication in each sampling point. That is, letting J :=I : £2r — > C be 
the pseudo-polar Fourier transform of an N x N image / and w : Qr M+ be any 
suitable weight function on £2r, the values 

^.(0)1, 0)2) C02)-a/w(c0i7c(^ for an{(Oi, (02) e Qr 

need to be computed. 

Let us comment on why the square root of the weight is utilized. If the weights 
w satisfy the condition in Theorem[T| we obtain P*wP = Id, which can be written in 
a symmetric form as follows: (y/wP)* y/wP — Id. This form shows that the operator 
y/wP can be inverted by taking the adjoint (y^f )*. In other words, each image can 
be reconstructed from its weighted pseudo-polar Fourier transform by applying the 
adjoint of the weighted pseudo-polar Fourier transform. This issue will be discussed 
in further detail in Subsection l2.4.2l 



2.3 Digital Shearlets on Pseudo-Polar Grid 

We next aim at deriving a faithful digitization of the shearlet transform associated 
with a band-limited cone-adapted discrete shearlet system to the pseudo-polar grid. 
This would settle Step (S3). 

2.3.1 Preparation for Faithful Digitization 

For this, let us recall the definition of the discrete shearlet transform associated 
with (|2]); taking the particular form ([T]i of the shearlet ^ e L^iM?) into account. 



Digital Shearlet Transforms 

Restricting our attention to the cone ^£2\, we obtain 
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2m (a ^_jS kill,) 



/,2-^2v/i(4-^-^i)V/2(^ + 2^"|)Z^,,e' 

for scale j, orientation k, position m, and cone I. 

To approach a faithful digitization, we first have to partition according to the 
partitioning of the plane into 'io\2, "^21, and "^22, as well as a centered rectangle 
31. The center as defined in (|6| will play the role of 31. Thus it remains to partition 
the set Q.R beyond the already defined partitioning into Q.\ and d\ (cf. Q and (|5]l) 
by setting 



12^ =X2^' U'^UXa]^ and ^2^ = ^2^' U "^U 12^- 



21, 



,22 



where 





= {(-1 


21; 

N 


2n \ 
R 1 


N 
2 


<i< 


N 1 ^ ,^ ^ RN\ 






= {i-T 


2e 

N 


2n \ 
R I 


N 
2 


<l< 


2 ' 2 — " — 


1}, 






2n 
R 


2^) 
N ) 


N 
2 


<l< 


N 1 ^ „ ^ RN\ 




^^R 




2n 
R 


2|) 
N ) 


N 
2 


<l< 


N _M <n<- 

2 ' 2 — " — 


!}■ 



When restricting to the cone i2|', say, the exact digitization of the coefficients of 
the discrete shearlet system is 

^ 7(0)1 , 0}2)2-j^{SlA^-jC0)e-^'''(^4-jSk'",'^) 
^ 7(0)1, a)2)2-^'2W(4-^«,)y(fe + 2>^)e-2^'(Vj-^*'"''») 



((ai,(02)Gi3j' 



RN 
2 



£ £ JiWuC02)2-^^W{4-J^^)V{k-2J+^l)e- 



(13) 



t- 2 



where V and W as well as the ranges of j, k, and m are to be carefully chosen. 

Our main objective will be to achieve a digital shearlet transform, which is an 
isometry. This - as in the continuum domain situation - is equivalent to requiring 
the associated shearlet system to form a tight frame for functions 7 : — > C. For 
the convenience of the reader let us recall the notion of a Parseval frame in this 
particular situation. A sequence {(Px);l£A ~ being some indexing set - is a tight 
frame for all functions 7 : Qr — > C, if 



I 



£ 7(0)1, a)2)9A(«i,£%) 
(coi,co2)ei2fi 



L |7(0)i,0)2)p. 
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In the sequel we will define digital shearlets on X2|^ and extend the definition to 
the other cones by symmetry. 

2.3.2 Subband Windows on the Pseudo-Polar Grid 

We start by defining the scaUng function, which will depend on two functions Vq and 
Wq, and the generating digital shearlet, which will depend on again two functions V 
and W. Wo and W will be chosen to be Fourier transforms of wavelets, and Vq and 
V will be chosen to be 'bump' functions, paralleling the construction of classical 
shearlets. 

First, let Wq be the Fourier transform of the Meyer scaling function such that 

suppWbC [-1,1] and Wo(±l)=0, (14) 
and let Vq be a 'bump' function satisfying 

supp Vo C [-3/2, 3 /2] with Vb(^ ) = 1 for |^ | < 1 , § e M. 
Then we define the scaling function (p for the digital shearlet system to be 

For now, we define it in continuum domain, and will later restrict this function to 
the pseudo-polar grid. 

Let next W be the Fourier transform of the Meyer wavelet function satisfying the 
support constraints 

suppWC [-4,1/4] U [1/4,4] and W(±l/4) = W(±4) = 0, (15) 

as well as, choosing the lowest scale Jl to be ji:= — [log4(^/2)] , 

[1084^1 

|Wo(4"^^^)P+ £ |W(4-^'^)|2 = 1 forallj^j <iV, ^ eK. (16) 

j=jL 

We further choose V to be a 'bump' function satisfying 

supp yc [-1,1] and V(±1)=0, (17) 

as well as 

|y(<?-l)|2 + |y(<^)|2 + |y(<^ + l)|2^1 forall |<^| < 1, <^ eM. (18) 
Then the generating shearlet \ir for the digital shearlet system on i2| is defined as 

iK(^i,^2) = W(^im|), (^i,&)eM^ (19) 
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Notice that ([18) implies 

V 

^ \V{2i^^s)\^^l forall|^|<l,^eM;y>0, (20) 

s=-2i 

which will become important for the analysis of frame properties. For the particular 
choice of Vb, Wo, V, and W in .ShearLab, we refer to Subsection |2.3.7| 



2.3.3 Range of Parameters 

We from now on assume that R and are both positive even integers and that = 
2"o for some integer «o e N. This poses no restrictions, since both parameters can 
be enlarged to satisfy this condition. 

We start by analyzing the range of j. Recalling the definition of the shearlet i/a 
in ([T9| and the support properties of W and V in ( fTS) and ([TtJ, respectively, we 
observe that the digitized shearlet 

2-4 W (4-^' f )y (/fc - 2^'+ 1 ^ )e2;ri(m,s[A^_, fl,) (21 ) 

from ( |T3] l has radial support 

„=4-'-'f h=0,...,A^-^ -^f (22) 

on the cone f2|'. To determine the appropriate range of y, we will analyze the 
precise support in radial direction. If j < — [log(/?/2)], then « < 1, which corre- 
sponds to only one point - the origin - and is dealt with by the scaling function. If 
j > [log4A^], we have « > Hence the value 1^(1/4) =0 (cf. (15\ ) is placed on 
the boundary, and these scales can be omitted. Therefore, the range of the scaling 
parameter will be chosen to be 

G {Jl, ■■■Jh}, where ;l := - [log(/?/2)] and jn := riog4A^] . 

Next, we determine the appropriate range of k. Again recalling the definition of 
the shearlet yf in ^T9\ , the digitized shearlet ( |2T| has angular support 

i^2-j-^N{k-l)+t2, t2=0,...,2-jN (23) 

on the cone i2|' . To compute the range of k, we start by examining the case j > 0. 
If k> 2j, we have £ > N/2. Hence the value V{-1) = (cf. ([T7]i) is placed on 
the seam line, and these parameters can be omitted. By symmetry, we also obtain 
k > —2^. Thus the shearing parameter will be chosen to be 



^e{-2^...,2^}. 
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2.3.4 Support Size of Shearlets 



We next compute the support as well as the support size of scaled and sheared ver- 
sion of digital shearlets. This will be used for the normalization of digital shearlets. 

As before, we first analyze the radial support. By ( [22] ), the radial supports of the 
windows associated with scales ji < j < Jh is 



n = 4^-'f+fi, fi=0, 



Aj-l . 15R 
I ^ 2 ' 



and the radial support of the windows associated with the scale Jl 
and Jh = riog4A^l are 



fi = 1, 
+ fi, ti=0, 



4^+1 R 

,f-4^«-lf, for;-;„. 



(24) 

Tlog4(^/2)l 
(25) 



Turning to the angular direction, by ( [23) , the angular support of the windows at 
scale j associated with shears — 2-' <k <2^ is 

i = 2--'-^N{k-i)+t2, t2=0,...,2-^'N, (26) 

the angular support at scale j associated with the shear parameter k = —2 ' is 

£ = 2-j-^N{-2j - 1) +f2, f2 = 2-^'f , . . .,2-jN, 

and for k — 2 ' it is 

i^2-J-^N{2j^l)+t2, f2=0,...,2-^f. (27) 

For the case < 0, we simply let ^ = and £ — —N/2 + 12 with t2 — 0,...,N. Also, 
for this lower frequency case, the window function W (4^-' (Bi)y(^ + 2-'^) is slightly 
modified to be W{4-j (Oi )Vo{k + 2j^). 

These computations now allow us to determine the support size of the function 
W (4->t0i )y {k + 2-1^) in terms of pairs {n,£), which for scale j and shear k, is 



^1 = 



4;-i . m 



1 



Jl <J< Jh, 
j = Jh, 



and 



2-^A^+l 



c^2 



2-q- 

N+l 



1 



-V <k< V with ./■ > 0, 
fee {-2^2'} with ;>0, 
i<0. 



(28) 



(29) 
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2.3.5 Digitization of the Exponential Term 

We next digitize the exponential term in ( |2T| i, which can be rewritten as 
We observe two obstacles: 

• The change of variables T := SlA^i^-jd) possible in ( [TJI can not be performed 
similarly in this situation due to the fact that the pseudo-polar grid is not invari- 
ant under the action of S^A^-j . This is however the first step in the continuum 
domain reasoning for tightness; see Chapter 1 1]. 

• The Fourier transform of a function defined on the pseudo-polar grid does not 
satisfy any Plancherel theorem. 

These problems require a slight adjustment of the exponential term, which will be 
the only adaption we allow us to make when digitizing. This will circumvent the 
two obstacles and enable us to construct a Parseval frame as well as derive a direct 
application of the inverse Fast Fourier transform in FDST 

The adjustment will be made by using the mapping : M \ {0} — > M defined by 
0{x^y) — (x, ^-). This yields the modified exponential term 

^-2;r,(m.(eo(5[)-l)(4-^-2^.4-^^t-2-^^)) ^ ^-2;r,(m.(4-^- ^,-2^+1 ^)) ^ 

which can be rewritten as 

^-2;r,-(m,(4-^-^,-2^+l^)) ^^-2;r,-(^+(l-^)„,2)^-2^M.(4-^-^,-2^+':i)>^ 

with t\ and f2 ranging over an appropriate set defined by ( |24] i, ( [ZS] ), and (|26ll-(p7]). 
Fig. |3]illustrates this adjustment. 



Fig. 3 Adjustment of the exponential tenn through the map 6 o {5j ) 



Now, taking into account of the support size of each W{A ■'Oi)y(A: + 2'^) as 
given in ( |28] l and ( p9| l, we obtain the following reformulation of ([30|: 

exp|-27ri( m, I -^-^^T fi, — ^^-^5- ^^lyj, fi,f2- (31) 
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This version shows that we might regard the exponential terms as characters of a 
suitable locally compact abelian group (see ifTSll ) with associated annihilator identi- 
fied with the rectangle 




where if/ and ^Jj^ were defined in ( |28] l and p9l ), respectively. This viewpoint 
will be crucial to guarantee that the digital shearlet system defined in Subsection 
2.3.6 provides a Parseval frame on the pseudo-polar grid Q^. In practice, pT| ) also 



ensures that in Step (S3) on each windowed image on the pseudo-polar grid only 
a 2D-iFFT - in contrast to a fractional Fourier transform - needs to be performed, 
thereby reducing the computational complexity. 

For the low frequency square, we further require the set 

^ = {(n,r2):n=-l,...,l,r2 = -f,...,f}, 

which will be shown to be sufficient for guaranteeing that digital shearlet system 
forms a Parseval frame. 



2.3.6 Digital Shearlets 

We are now ready to define digital shearlets, which we define as functions on the 
pseudo-polar grid Qr. The spatial domain picture can thus be derived by the inverse 
pseudo-polar Fourier transform. 



Definition 2. Retaining the definitions and notations from Subsection |2.3| for all 

(01,02) e ^ff', we define digital shearlets at scale j e {^'l, • • • ,7//}, shear k = 
[—2-' , 2-'] n Z, and spatial position m G by 

CTli.J«i,02) = ^H^(4->Oi)W(fc + 2>^)X^.(oi,02)/K-^^^^^ 
where Y' = V for ; > and Y' = Yo for ; < 0, and 



C(0i,02) = 




(01,02) ^^s'u^i, 
(oi,02)e (^K'u^i)\'^, 
(01,02) e 'if. 



The shearlets crjj, C7/|,„ on the remaining cones are defined accordingly by 

symmetry with equal indexing sets for scale j, shear k, and spatial location m. For 
Iq — 1,2, (oi , O2) G £2^, and no G we define the scaling function 

<(0i,02) = ^M^(cOi,fi)2)Z^.o(Oi,02)^.2^«o,(?.^)). 
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Then the digital shearlet system DSH is defined by 

DSH = {(p,;« : lo = l,2,«o e ^}U {cjj_,^„, : ; € {jt, . . .,;//}, ^ e {-2\2^}, 

me^j^„i = 11,12,21,22}. 

As desired, the digital shearlet system DSH, which we derived as a faithful digiti- 
zation of the continuum domain band-limited cone-adapted discrete shearlet system, 
forms a Parseval frame for J : Qr — > C. 

Theorem 2 (f201). The digital shearlet system DSH defined in Definition^^orms a 
Parseval frame for functions J : Qg C. 

Proof. Letting J : £2r C, we claim that 



k)-'^0 lj,k,m 



(32) 



which proves the result. Here {Jy^JijaR L(<Bi,ffl2)Gi2«-/i(«i,fi)2)-/2(«i,ft>2) for 
71,72 :i2R^C. 

We start by analyzing the first term on the RHS of ( (32] i. Let lo € {1,2} and 
7c : C be defined hy J c{o}\,0}2) := C{o}i,0}2) ■ J {o}\,o>i) for (coi , 0)2) e 12^. 

Using the support conditions of ^, 



«0 



7(t(Ji, 0)2)^^0 (0)1, 0)2) 



= WX I E /c(«l,tO2)-0(tOi,tO2)-.-2-(«O.(f.^)) 

I-^I no «=-l(!=-7V/2 

The choice of now allows us to use the Plancherel formula, see Subsection [533] 
Exploiting again support properties (see Subsection 2.3.5| l, we conclude that 

Y.\^J,<')n,?^ E |C(a)i,a)2)-7(coi,co2)|2.|(^(«i,co2)|2. 



Combining Iq = 1,2 and using ( [T4] i, we proved 



EEl(-/,Oi3j 

lo "0 



^ 17(0)1, a)2)|2.|Wo(a)i)|2. 

((Oi,6l>2)Gi2/i 



(33) 



Next we study the second term on the RHS in ([32J. By symmetry, it suffices to 
consider the case i =21. By the support conditions on W and V (see ([TSj and (fTTjl), 



E K-^'^i.L")^J^^E E E ■/(«1,»2)ct2[,„(C0i,0)2; 
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1 



= Lt^i L L Jc{cOhO>2)-Wi4-jcoi) 
j,k \'^jMmeMj^t {(Ouc>i)enf 

^ 4J+'(«/2) 2-i-^N{k+l) 

= T^w^\ EL L -/c(«i,o>2) 

},k \'^''JM me.'^j^t n=4J-i{R/2)e=2-J-'N{k-l) 

■W{4-J(0i) ■ Vj{k + 2.>^) ■ 6-2^^' ('-(4-^' 1? .-2^'+' ^ )> 

Similarly as before, the choice of ^j^k does allow us to use the Plancherel formula, 
see Subsection |2.3.5| Hence, 



Next we use (|20ll to obtain 



Jc{(Oi,0)2)-W{4-jcoi)Vj{k + 2j 



CD, 



E Jc{couO>2)-W{4-Jcoi)-VJ{k + 2J^) 



((B,,(02)Gi2^' 



.7=;l 

JH 



k=-V 



((Bl,(02)Gi2|' i=}L 

Thus the second term on the RHS in (|32]) equals 



I L l(-/,cTj,,,J^,|2= L \J{(oum)?- L .|w(4-^«i)|2 

Finally, our claim ( (32] l follows from combining ( |33| l, p4l ), and ( [T6] l. □ 



(34) 



2.3.7 Digital Shearlet Windowing 

The final Step (S3) of the FDST then consists in decomposing the data on the points 
of the pseudo-polar grid given by the previously - in Steps (SI) and (S2) - computed 
weighted pseudo-polar image J^. : Qr — > C into rectangular subband windows ac- 
cording to the digital shearlet system DSH defined in Definition |2] followed by a 
2D-iFFT More precisely, given J^,, the set of digital shearlet coefficients 



4o:=(^,<o)c2„ forallio,«o 
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and 



'j.k.m 



j.k,m 



for all j,k,m, i 



is computed followed by application of the 2D-iFFT to each windowed image J„, (Pq 
and -/wCj restricted on the support of (Pq° and aj ^ q, respectively. 

The definition of the digital shearlet system DSH in Definition|2]requires appro- 
priate choices of the functions 0, Vq, V, Wq, and W, and the required conditions are 



stated throughout Subsection 2.3.2 We now discuss one particular choice, which is 
chosen in ShearLab We start selecting the 'wavelets' Wo and W. In Subsection 



2.3.2 these functions were defined to be Fourier transforms of the Meyer scaling 
function and wavelet function, respectively, i.e.. 



1 



Wo(^) = <! cos[fv(i|^|-i)] 




\^\<l 

l<\^\<h 

otherwise, 



and 



sin[fv(f|^|-i)] 
W(^)=<; cos[fv(i|^|-i)] 




3< 



1<I^I<4, 
otherwise, 



where v > is a C*^ function or C°° function such that v{x) + v{l — x) = 1 for < 
x<l. One possible choice for v is the function v{x) = x^{35 — S4x + 70x^ — 2Qx^), 
< X < 1, which then automatically fixes Wo and W. Since |Wo(<^ ) p + |W(^ ) ^ = 1 
for \ < 1, the required condition ([16]) is satisfied. The graphs of this choice of 
functions Wo, W, and v are illustrated in Fig. [4] 





Wo 



W 



Fig. 4 The graphs of v, Wq, and w. 



The function v can be also used to design the 'bump' function V as well, 
which needs to satisfy ( [T8] l. One possible choice for V is to define it by V{£,) — 
^/v{Y+~^Y+V(T^^^, — 1 < < 1. Vb can then simply be chosen as Vo = 1- 

Let us finally mention that (j) is defined depending on Vq and Wo, wherefore fixing 
these two functions determines (j) uniquely. 
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2.4 Algorithmic Realization of the FDST 



We have previously discussed all main ingredients of the fast digital shearlet trans- 
form (FDST) - Fast PPFT, Weighting, and Digital Shearlet Windowing -, and will 
now summarize those findings. Depending on the application at hand, a fast inverse 
transform is required, which we will also detail in the sequel. In fact, we will present 
two possibilities; the Adjoint FDST and the Inverse FDST depending on whether the 
weighting allows to use the adjoint for reconstruction or whether an iterative proce- 
dure is required for higher accuracy. Fig. |5]provides an overview of the main steps 
of of the FDST and its inverse. For a more detailed description of FDST, Adjoint 
FDST, and Inverse FDST in form of pseudo-code, we refer to I.20J . 




Fig. 5 Flowcharts of the FDST (left) and its inverse (right). 



For the sake of brevity, we now let P, w, and W denote the Fast PPFT from 
Subsection 2.1.2 the weighting on the pseudo-polar grid described in Subsection 



2.2.4 and windowing operator consisting of the application of the shearlet windows 



followed by 2D-iFFT to each array as detailed in Subsection |2. 3. 7[ respectively. 



2.4.1 FDST 



We can summarize the steps of the algorithm FDST as follows: 

• Step (SI): For a given image /, apply the Fast PPFT as described in Subsection 
|2.1.2| to obtain the function PI :£1r~^ C. 

• Step (S2): Apply the square root of an off-line computed weight function w : 



Q.R C to PI as described in Subsection 2.2.4 yielding ^JwPl : Qr C. 
Step (S3): Apply the shearlet windows to the function wPI, followed by a 2D- 
iFFT to each array to obtain the shearlet coefficients Wy/wPI, which we denote 

by c?(,, io,no and c^^, ,,,, j,k,m,i. 



Digital Shearlet Transforms 



25 



2,4.2 Adjoint FDST 

Assuming that the weight function w used in Step (S2) satisfies the condition in 
Theorem [T] and using the Parseval frame property of the digital shearlet system 
(Theorem |2]i, we obtain 

{Wy/^pyw^wP = p*^{w*w)^p = p'-wP = Id. 

Hence in this case, the FDST, which is abbreviated by Wy/wP can be inverted by 
applying the adjoint FDST, which cascades the following steps: 

• Step 1: For given shearlet coefficients C, i.e., c,'", lo,no and cjj^,„, j,k,m,i, 
compute the linear combination of the shearlet windows with coefficients cj,", 
lo,no and c^-^^^^, j,k,m,i. This gives the function W*C : Qr — C. 

• Step 2: Apply the square root of an off-line computed weight function w : Qr — > 
C to W*C, yielding the function y^W*C -.Qr^C. 

• Step 3: Apply the Fast Adjoint PPFT by running the Fast PPFT 'backwards'. 
For this, we just notice that the adjoint fractional Fourier transform of a vector 
c G C'^+' with respect to a constant a e C is given by F^'^^c. Also, for m > 
N, the adjoint padding operator ^ applied to a vector c e C™ is given by 
{E*^ f^c)(k) = c{k), k = -N/2, . . . ,N/2 - 1. The Adjoint PPFT gives an image 



2.4.3 Inverse FDST 



Normally - as also with the relaxed form of weights debated in Subsection 2.2.2 
the weights will not satisfy the conditions of Theorem [T] precisely. A measure for 
whether application of the adjoint is still feasible was akeady discussed in Subsec- 
tion |2]23] (see also Subsection |4.2| i. If higher accuracy of the reconstruction is re- 
quired, one might use iterative methods, such as conjugate gradient methods. Since 
the digital shearlet system forms a Parseval frame, we always have 

w-wVn'P = VwP 



Hence, iterative methods need to be 'only' applied to reconstruct an image / from 
knowledge of J := y/wPI, i.e., to solve the equation 

P*wPI = P*wJ 

for /. Since J might not be in the range of P, I is typically computed by solving 
the weighted least square problem min/ 1| y^P/ — y^7||2. Since the matrix cor- 
responding to P*P is symmetric positive definite, iterative methods such as the 
conjugate gradient method are applicable. The conjugate gradient method is then 
applied to the equation Ax = b with A = P*wP and b — P*wJ. Its performance 
can be measured by the condition number of the operator P*wP: cond{P*wP) — 
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^nmx{P*wP) / X,ni„{P*wP), and it turns out that the weight function serves as a pre- 



conditioner. We remark that this measure is more closely studied in Subsection 4.2 

To illustrate the behavior of the weights with respect to this measure, in Table]?] 
we compute cond{P*wP) for the weight functions arising from Choices 1 and 2 (cf. 
Subsection 2.2.2| i with oversampling rate = 8. Notice that the condition numbers 
of P*wP are generally smaller than 2. 



Table 3 Comparison of Choices 1 and 2 based on the performance measure cond{P*wP). 



N 


32 


64 


128 


256 


512 


Choice 1 


1.379 


1.503 


1.621 


1.731 


1.833 


Choice 2 


1.760 


1.887 


2.001 


2.104 


N/A 



3 Digital Shearlet Transform using Compactly Supported 
Shearlets 

In this section, we will discuss two implementation strategies for computing shearlet 
coefficients associated with a cone-adapted discrete shearlet system now based on 
compactly supported shearlets, as introduced in Chapter 1 1 1. Again, one main focus 
will be on deriving a digitization which is faithful to the continuum setting. 

Recall that in the context of wavelet theory, faithful digitization is achieved by 
the concept of multiresolution analysis, where scaling and translation are digitized 
by discrete operations: Downsampling, upsampling and convolution. In the case of 
directional transforms however, three types of operators: Scaling, translation and 
direction, need to be digitized. In this section, we will pay particular attention to 
deriving a framework in which each of the three operators is faithfully interpreted 
as a digitized operation in digital domain. Both approaches will be based on the 
following digitization strategies: 

• Scaling and translation: A multiresolution analysis associated with anisotropic 
scaling can be applied for each shear parameter k. 

• Directionality: A faithful digitization of shear operator -Sj ;/!^ has to be achieved 
with particular care. 

After stating and discussing the two main obstacles we are facing when considering 



compactly supported shearlets in Subsection 3.1 we present the digital separable 
shearlet transform (DSST), which is associated with a shearlet system generated 
by a separable function alongside with discussions on its properties, e.g., its redun- 



dancy; see Subsection 3.2 Subsection 3.3 then presents the digital non-separable 



shearlet transform (DNST), whose shearlet elements are generated by non- separa- 
ble shearlet generator. 
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3.1 Problems with Digitization of Compactly Supported Shearlets 



Compactly supported shearlets have several advantages, and we exemplarily men- 
tion superior spatial localization and simplified boundary adaptation. However, we 
have to face the following two problems: 

(PI) Compactly supported shearlets do not form a tight frame, which prevents uti- 
lization of the adjoint as inverse transform. 

(P2) There does not exist a natural hierarchical structure, mainly due to the appli- 
cation of a shear matrix, which - unlike for the wavelet transform - does not 
allow a multiresolution analysis without destroying a faithful adaption of the 
continuum setting. 

Let us now comment on these two obstacles, before delving into the details of 



the implementation in Subsection 3.2 



3.1.1 Tightness 

Let us first comment on the problem of non-tightness. Letting ((7,),g/ denote a frame 
for L^(R^) - for example, a shearlet frame -, each function / e L^(M^) can be 
reconstructed from its frame coefficients ((/, CJ,)),^/ by 

f = adS-\a.), 
iei 

where S — T^iei {' ^ (^i) (^i is the associated frame operator on L^(M^), see Chapter (Tl. 
However, in case that ((T,),g/ does not form a tight frame, it is in general difficult to 
explicitly compute the dual frame elements (ff,). 

Nevertheless, it is well known that the inverse frame operator S^^ can be effec- 
tively approximated using iterative schemes such as the Conjugate Gradient method 
provided that the frame (c7,),g/ has 'good' frame bounds in the sense of their ratio 
being 'close' to 1, see also (2T\. Therefore, now focussing on the situation of shear- 
let frames, we may argue that input data / can be efficiently reconstructed from its 
shearlet coefficients, if we use a compactly supported shearlet frame with 'good' 
frame bounds. In fact, the theoretical frame bounds of compactly supported shearlet 
frames have been theoretically estimated as well as numerically computed in L19l . 
These results were derived for the class of 2D separable shearlet generators y/ al- 
ready described in Chapter (T\, which we briefly recall for the convenience of the 
reader: 

For positive integers K and L fixed, let the ID lowpass filter mo be defined by 

|mo(^i)|2 = (cos(7r^,))'^ L ^ (M^^i)?", 
« =0 V n / 

for e M. Further, define the associated bandpass filter mi by 
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|mi(^l)l'-K(^i + l/2)P, ^leM, 
and the ID scaling function by 

1=0 

Using the fiher mi and scaHng function 0i, we now define the 2D scaling function 
^ and separable shearlet generator xjf by 

= and v>(^i,fe)=mi(4<^i)0i(^i)0,(2&). (35) 

In |fT9l . it was shown that compactly supported shearlets yfj^k.m generated by the 
shearlet generator form a frame for L^(C)^ with appropriately chosen parameters 
K and L, where 

C = {^eR^:\^2/^i\<l, |^i|>l}. 

This construction is directly extended to construct a cone-adapted discrete shearlet 
frame fov L^{R^) (cf. also Chapters [1 1 and |4|). 

Table |4]provides some numerically estimated frame bounds in L^{C) for certain 
choice of K and L. It shows that indeed the ratio of the frame bounds of this class 
of compactly supported shearlet frames is sufficient small for utilizing an iterative 
scheme for efficient reconstruction; in this sense the frame bounds are 'good'. 



Table 4 Numerically estimated frame bounds for various choices of the parameters K and L. ci 
and C2 are the sampling constants in the sampling matrix M^. for translation (see Chapter 1 1 1). 



K 


L 




C2 


B/A 


39 


19 


0.90 


0.15 


4.1084 


39 


19 


0.90 


0.20 


4.1085 


39 


19 


0.90 


0.25 


4.1104 


39 


19 


0.90 


0.30 


4.1328 


39 


19 


0.90 


0.40 


5.2495 



The frequency covering by compactly supported shearlets \j/j,k.m7 
+ L E |v>(5[A2,^)|2 + |f(5[A2;^)|2, 

j>OkeKj 

is closely related to the ratio of frames bounds and, in particular, which areas in 
frequency domain cause a larger ratio. This function is illustrated in Fig. |6j which 
shows that its upper and lower bounds are as expected well controlled. 
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(a) Whole frequency plane. (b) Horizontal cone. (c) Vertical cone. 

Fig. 6 Frequency covering by shearlets |v/y,j:.„,p: (a) Frequency covering of the entire frequency 
plane, (b) Frequency covering of the horizontal cone, (c) Frequency covering of the vertical cone. 

3.1.2 Hierarchical Structure 

Let us finally comment on the problem to achieve a hierarchical structuring. To al- 
low fast implementations, the data structure of the transform is essential. The hierar- 
chical structure of the wavelet transform associated with a multiresolution analysis, 
for instance, enables a fast implementation based on filterbanks. In addition, such 
a hierarchical ordering provides a full tree structure across scales, which is of par- 
ticular importance for various applications such as image compression and adaptive 
PDE schemes. It is in fact mainly due to this property - and the unified treatment of 
the continuum and digital setting - that the wavelet transform became an extremely 
successful methodology for many practical applications. 

From a certain viewpoint, shearlets i/^.A ,m can essentially be regarded as wavelets 
associated with an anisotropic scale matrix , when the shear parameter k is fixed. 
This observation allows to apply the wavelet transform to compute the shearlet co- 
efficients, once the shear operation is computed for each shear parameter k. This 
approach will be undertaken in the digital formulation of the compactly supported 
shearlet transform, and, in fact, this approach implements a hierarchical structure 
into the shearlet transform. The reader should note that this approach does not lead 
to a completely hierarchical structured shearlet transform - also compare our dis- 
cussion at the beginning of this section -, but it will be sufficient for deriving a fast 
implementation while retaining a faithful digitization. 



3.2 Digital Separable Shearlet Transform (DSST) 

We now describe a faithful digitization of the continuum domain shearlet transform 
based on compactly supported shearlets as introduced in [22J, which moreover is 
highly computationally efficient. 
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3.2.1 Faithful Digitization of the Compactly Supported Shearlet Transform 

We start by discussing those theoretical aspects which allow a faithful digitization 
of the shearlet transform associated with the shearlet system generated by ( (35] l. For 
this, we will only consider shearlets \l/j^k.m for the horizontal cone, i.e., belonging to 
f (i/,c). Notice that the same procedure can be applied to compute the shearlet co- 
efficients for the vertical cone, i.e., those belonging to 'f'( V/, c), except for switching 
the order of variables. 

To construct a separable shearlet generator y/ G L^(M^) and an associated scaling 
function (j) G L^(M^), let G L^(M) be a compactly supported ID scaling function 
satisfying 

(l)i{xi) = h{ni)V2^i{2xi-ni) (36) 

for some 'appropriately chosen' filter h - we comment on the required condition be- 
low. An associated compactly supported ID wavelet Xj/i G L^(E) can then be defined 
by 

Vi(^i)= E «("i)V20i(2xi-ni), (37) 

n,eZ 

where again g is an 'appropriately chosen' filter The selected shearlet generator is 
then defined to be 

V/(xi,X2) = V^i(jci)0i(x2), (38) 

and the scaling function by 

(j){xi,X2) = 0i(jCi)0i(x2). 

Let us comment on whether this is indeed a special case of the shearlet generators 
defined in ( |35| ). The Fourier transform of y/ defined in ( |38| ) takes the form 

V/(^i,^2)=mi(^i/2)0i(^i/2)0i(^2/2), 

where mi is a trigonometric polynomial whose Fourier coefficients are g{ni). We 
need to compare this expression with the Fourier transform of the shearlet generator 
\j/ given in ([35|, which is 

VA(^,,^2)=mi(4^i)0i(25)^i(^2), 

with ID scaling function 01 defined in ( (36] l. We remark that this later scaling func- 
tion is slightly different defined as in ( |35| ). This small adaption is for the sake of 
presenting a simpler version of the implementation; essentially the same implemen- 
tation strategy as the one we will describe can be applied to the shearlet generator 
given in ( |35] l. 

The filter coefficients h and g are required to be chosen so that y/ satisfies a 
certain decay condition (cf lfT9l of Chapter LU) to guarantee a stable reconstruction 
from the shearlet coefficients. 
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For the signal / e L^(M^) to be analyzed, we now assume that, for 7 > fixed, / 
is of the form 

/W = L fjin)2'^i2-'xi-nu2'x2-n2). (39) 

neZ^ 

Let us mention that this is a very natural assumption for a digital implementation in 
the sense that the scaling coefficients can be viewed as sample values of / - in fact 
fj{n) = /(2^"'n) with appropriately chosen (j). Now aiming towards a faithful digi- 
tization of the shearlet coefficients (/, y/j,ic,m) for j — 0,. . . ,J ~ I, we first observe 
that 

(/, Wj.k.,n) = (/(52-,/2,(-)), </;,().„,(•)), (40) 

and, WLOG we will from now on assume that j/2 is integer; otherwise either \j/2'] 
or [j/2\ would need to be taken. Our observation ( |40l ) shows us in fact precisely 
how to digitize the shearlet coefficients (/. yfj,k,m)- By applying the discrete sepa- 
rable wavelet transform associated with the anisotropic sampling matrix to the 
sheared version of the data f {82-] /i /,{■))■ This however requires - compare the as- 
sumed form of / given in ( |39] l - that /(5'2-j/2j-(-)) is contained in the scaling space 

Vj = {2-'<l) (2-' • ~«i , 2-' • -«2) : («i , «2) e 

It is easy to see that, for instance, if the shear parameter 2^^l^k is non-integer, this 
is unfortunately not the case. The true reason for this failure is that the shear matrix 
52-j/2^ does not preserve the regular grid 2^''l? in Vj, i.e.. 

In order to resolve this issue, we consider the new scaling space V'y^_y/2 j defined by 

yljIU = {2^+^/''0(5,(2^+^-/2 . . _„2)) : («i ,«2) G I?}. 

We remark that the scaUng space yj^ji2 j is obtained by refining the regular grid 
2~''l? along the jcpaxis by a factor of 2^1^. With this modification, the new grid 
2^''^-'/^Z X 2^''Z is now invariant under the shear operator 5'2-j/2^, since with Q = 
diag(2,l), 

2-^-^l^'L X 2-^1 = 2-^Q-'l^{I?) ^ 2-^Q-'l^{Sk{l})) 
= 52-;/2^(2--'->/2zx2--'Z). 

This allows us to rewrite f{S2-j/2i^{-)) in ( |40j i in the following way. 

Lemma 1. Retaining the notations and definitions from this subsection, letting f 
2'/^ and *i denote the ID upsampling operator by a factor 0/ 2'/^ and the ID 
convolution operator along the x\-axis, respectively, and setting hji2{ni) to be the 
Fourier coefficients of the trigonometric polynomial 
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%2(<^i)='n' E Kn,)e-^-'^'"^^\ (41) 

we obtain 

f{S2-m{^))^ L fj{Skn)2'+^'^h{^'+'''^x,-n,,2'x2-n2), 

nel? 

where 

fj{n) = ((//)|2J/2*i^y2)(«)- 

The proof of this lemma requires the following result, which follows from the 
cascade algorithm in the theory of wavelet. 

Proposition 1 ( II22I ). Assume that (j)i and Xj/i e L^(M) satisfy equations ^36\ and 
( |37| l respectively. For positive integers jy < j2, we then have 

2T0i(2^ixi-«i)= ^ (fl'i-2^2->ini)2 * 01 (2>2xi-rfi) (42) 

d^ez 

and 

2Tv/i(2^ixi-«i)= 8j,-j,idi-2j^--j'ni)2^^i{2j^xi-di), (43) 

diez 

where hj and gj are the Fourier coefficients of the trigonometric polynomials Hj 
defined in ( |41| l and Gj defined by 

k=0niez ineZ 

for i > Q fixed. 

Proof (Proof of Lemma^. Equation ( |42j i with j\ = J and j2 =J + j /2 implies that 

2^/201(2^X1 -m) = V,./2(^/i -2^/2„i)2^/2+VVi(2'+^/'^i -^O- (44) 
dieZ 

Also, since is a 2D separable function of the form 0(xi,X2) = 0i (xi)0i (^2), we 
have that 

/W= E (E /^(«i'«2)2^/'0i(2-'xi-«i))2^/20i(2^X2-«2). 
By ( [44l (, we obtain 

.m= E /y(«)2'+^^/^0(2^e^^/2x-n), 
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where Q — diag(2, 1). Using Q^I^S2-ji2i^ — S^Q^^^, this finally implies 

nel? 

= £ /,(«)2^+^-/4^(5^(2/g;72^_5_,„)) 

nel? 

The lemma is proved. □ 

The second term to be digitized in ( |40] i is the shearlet VO' ^-.m itself. A direct corol- 
lary from Proposition [T| is the following result. 

Lemma 2. Retaining the notations and definitions from this subsection, we obtain 
WiXmix) = E gJ-Mi-2'-^m,)hj_^,2{d2^2'-'l^m2)2'+'l^<p{2'Qil^x-d). 

As already indicated before, we will make use of the discrete separable wavelet 
transform associated with an anisotropic scaling matrix, which, for ji and j2 > as 
well as c e ^(Z^), we define by 

^ji..h[c){ni,n2)^ E gyj(mi -2-"ni)/i^2(m2-2''2«2)c(mi,m2), (ni,n2)eZ^- 

(45) 

Finally, Lemmata [T] and [2] yield the following digitizable form of the shearlet 
coefficients {f,y/j^k,m)- 

Theorem 3 ( II22I ). Retaining the notations and definitions from this subsection, and 
letting], 2'/^ be ID downsampling by a factor of 2^1^ along the horizontal axis, we 
obtain 

if, WjX'n) - Wj_jj_jf2[(^{fj{Sk-) * *1 V)^2Jy2) 

where <i>k{n) = {^{Sk{-)),^{- - n)) for n e I?, andhji2{n\)^hji2{-n\). 
3.2.2 Algorithmic ReaUzation 

Computing the shearlet coefficients using Theorem [3] now restricts to applying the 
discrete separable wavelet transform ( |45] l associated with the sampling matrix A2j 
to the scaling coefficients 

5^-;/2, (//)(«) ((/y(5r) * *i ^/2)^2^.^,(«) for fj e (46) 

Before we state the explicit steps necessary to achieve this, let us take a closer look 
at the scaling coefficients S'^_j^2i^{fj)' which can be regarded as a new sampling of 
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Fig. 7 Illustration of appli- 
cation of the digital shear 
operator S'^^^^: The dashed 
lines correspond to the refine- 
ment of the integer grid. The 
new sample values lie on the 
intersections of the sheared 
lines associated with Si /4 with 
this refined grid. 



the data fj on the integer grid I? by the digital shear operator S'^2-Jftk' procedure 
is illustrated in Fig.jvjin the case l^^l^k = — 1/4. 

Let us also mention that the filter coefficients ^k{n) in (|46| can in fact be easily 
precomputed for each shear parameter k. For a practical implementation, one may 
sometimes even skip this additional convolution step assuming that = X{Qia)- 

Concluding, the implementation strategy for the DSST cascades the following 
steps: 

• Step 1: For given input data fj, apply the ID upsampling operator by a factor 
of 2^1^ at the finest scale j = J. 

• Step 2: Apply ID convolution to the upsampled input data fj with ID lowpass 
filter hji2 at the finest scale /' = J. This gives fj. 

• Step 3: Resample fj to obtain fj{Sk{n)) according to the shear sampling matrix 
Sk at the finest scale j = J. Note that this resampling step is straightforward, 
since the integer grid is invariant under the shear matrix Sk- 

• Step 4: Apply ID convolution to fj{Sk{n)) with hjji followed by ID downsam- 
pling by a factor of 2-'/^ at the finest scale j — J. 

• Step 5: Apply the separable wavelet transform W/_y7_j/2 across scales j = 
0,1,. ..,7-1. 

3.2.3 Digital Realization of Directionality 

Since the digital realization of a shear matrix 52-j/2^ by the digital shear operator 
si /, , is crucial for deriving a faithful digitization of the continuum domain shearlet 
transform, we will devote this subsection to a closer analysis. 

We start by remarking that in fact in the continuum domain, at least two opera- 
tors exist which naturally provide directionality: Rotation and shearing. Rotation is 
a very convenient tool to provide directionality in the sense that it preserves impor- 
tant geometric information such as length, angles, and parallelism. However, this 
operator does not preserve the integer lattice, which causes severe problems for dig- 
itization. In contrast to this, a shear matrix Sk does not only provide directionality, 
but also preserves the integer lattice when the shear parameter k is integer Thus, it 
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is conceivable to assume that directionality can be naturally discretized by using a 
shear matrix S^- 

To start our analysis of the relation between a shear matrix S2-j/2i^ and the asso- 
ciated digital shear operator S'^-j/Tf.^ let us consider the following simple example: 
Set fc = X{x:xt=o}- Then digitize fc to obtain a function defined on I? by set- 
ting fd{n) = fc{n) for all n € Z^. For fixed shear parameter i G M, apply the shear 
transform Ss to fc yielding the sheared function fc{Ss{ )). Next, digitize also this 
function by considering fc{Ss{-))\i2- The functions and fc{Ss{-))\i2 are illus- 
trated in Fig. |8]for 5 = — 1 /4. We now focus on the problem that the integer lattice is 




not invariant under the shear matrix 5'i/4. This prevents the sampling points S\n{n), 
n e from lying on the integer grid, which causes aliasing of the digitized im- 
age fc{S_\i/^{-))\j2 as illustrated in Fig.lola). In order to avoid this aliasing effect, 
the grid needs to be refined by a factor of 4 along the horizontal axis followed by 
computing sample values on this refined grid. 

More generally, when the shear parameter is given one can 

essentially avoid this directional aliasing effect by refining a grid by a factor of 
2^/^ along the horizontal axis followed by computing interpolated sample values on 
this refined grid. This ensures that the resulting grid contains the sampling points 
((2-^/2^)«2,«2) for any «2 G ^ and is preserved by the shear matrix S_2-j/2i^- 
This procedure precisely coincides with the application of the digita l shear oper- 

in which the 



ator S'^-jfti^y i-S-, we just described Steps 1-4 from Subsection 3.2.2 
new scaling coefficients S'^_jp^{ fj){n) are computed. 

Let us come back to the exemplary situation of fc = X{x:xi=o} and 5_i/4 we 
started our excursion with and compare /c('5-i/4('))lz2 with S'^^j^{fd)\x2 obtained 
by applying the digital shear operator S'^^^^ to And, in fact, the directional alias- 
ing effect on the digitized image fc{S-i/4{n)) in frequency illustrated in Fig.|9|a) 
is shown to be avoided in Fig.|9j(b)-(c) by considering S''_^^^{fii)\j2. Thus applica- 
tion of the digital shear operator S'^_jj2^, allows a faithful digitization of the shearing 
operator associated with the shear matrix S2-j/2i^- 
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3.2.4 Redundancy 

One of the main issues which practical appHcabihty requires is controllable redun- 
dancy. To quantify the redundancy of the discrete shearlet transform, we assume that 
the input data / is a finite linear combination of translates of a 2D scaling function 
([) at scale J as follows; 

2-'-i2-'-i 

« 1=0^2=0 

as it was already the hypothesis in ([39|. The redundancy - as we view it in our 
analysis - is then given by the number of shearlet elements necessary to represent 
/. Furthermore, to state the result in more generality, we allow an arbitrary sampling 
matrix = diag(ci ,C2) for translation, i.e., consider shearlet elements of the form 

VjAmi-) = 2^j\i/{SkA2j ■ -M,m). 

We then have the following result. 

Proposition 2 ( ll22l ). The redundancy of the DSST is 




Proof. For this, we first consider shearlet elements for the horizontal cone for a 
fixed scale j e {0, ... ,7 — 1}. We observe that there exist 2^/^+' shearing indices k 
and 2-' • 2-'/'^ • (ciC2)^' translation indices associated with the scaling matrix and 
the sampling matrix M^, respectively. Thus, 2^-'+' (ciC2)^' shearlet elements from 
the horizontal cone are required for representing /. Due to symmetry reasons, we 
requke the same number of shearlet elements from the vertical cone. Finally, about 
c^^ translates of the scaling function are necessary at the coarsest scale j = 0. 

Summarizing, the total number of necessary shearlet elements across all scales 
is about 
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The redundancy of each shearlet frame can now be computed as the ratio of the 
number of coefficients d„ and this number. Letting J ^ °o proves the claim. □ 

As an example, choose a translation grid with parameters ci = 1 and C2 = 0.4. 
Then the associated DSST has asymptotic redundancy 10/3. 



3.2.5 Computational Complexity 

A further essential characteristics is the computational complexity (see also Subsec- 
tion |4.6| l, which we now formally compute for the discrete shearlet transform. 

Proposition 3 ( II19I ). The computational complexity of the DSST is 

Q^2^og2il/2iL/2-l)) j^^j^y 

Proof. Analyzing Steps 1-5 from Subsection |3.2.2| we observe that the most time 
consuming step is the computation of the scaling coefficients in Steps 1-4 for the 
finest scale j — J. This step requires ID upsampling by a factor of 2^/^ followed 
by ID convolution for each direction associated with the shear parameter k. Letting 
L denote the total number of directions at the finest scale j = J, and the size of 
2D input data, the computational complexity for computing the scaling coefficients 
in Steps 1 - 4 is 0{2^/^L ■ N). The complexity of the discrete separable wavelet 
transform associated with for Step 5 requires 0{N) operations, wherefore it is 
negligible. The claim follows from the fact that L = 2{2-Vl^ + 1). □ 

It should be noted that the total computational cost depends on the number L of 
shear parameters at the finest scale j = J, and this total cost grows approximately by 
a factor of as L is increased. It should though be emphasized that L can be chosen 
in such a way that this shearlet transform is favorably comparable to other redun- 
dant directional transforms with respect to running time as well as performance. A 
reasonable number of directions at the finest scale is 6, in which case the constant 
factor 2'°S2(i/2(i-/2-i)) Propositionjsjequals 1. Hence in this case the running time 
of this shearlet transform is only about 6 times slower than the discrete orthogo- 
nal wavelet transform, thereby remains in the range of the running time of other 
directional transforms. 



3.2.6 Inverse DSST 



In Subsection 3. LI we already discussed that this transform is not an isometry, 
wherefore the adjoint cannot be used as an inverse transform. However, the 'good' 
ratio of the frame bounds in the sense as detailed in Subsection l3.1.1l leads to a fast 
convergence rate of iterative methods such as the conjugate gradient method. Let us 
mention that using the conjugate gradient method basically requires computing the 
forward DSST and its adjoint, and we refer to [24.1 and also Subsection 2.4.2 for 
more details. 
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3.3 Digital Non-Separable Shearlet Transform (DNST) 



In this section, we describe an alternative approach to derive a faithful digitaliza- 
tion of a discrete shearlet transform associated with compactly supported shearlets. 
This algorithmic realization, which was developed in ||23]| . resolves the following 
drawbacks of the DSST: 

• Since this transform is not based on a tight frame, an additional computational 
effort is necessary to approximate the inverse of the shearlet transform by itera- 
tive methods. 

• Computing the interpolated sampling values in ( |46l ) requires additional compu- 
tational costs. 

• This shearlet transform is not shift-variant, even when downsampling associated 
with is omitted. 

We emphasize that although this alternative approach resolves these problems, the 
algorithm DSST provides a much more faithful digitalization in the sense that the 
shearlet coefficients can be exactly computed in this framework. 

The main difference between DSST and DNST will be to exploit non-separable 
shearlet generators, which give more flexibility. 



3.3.1 Shearlet Generators 

We start by introducing the non-separable shearlet generators utilized in DNST. 
First, for each scale parameter j > 0, define the shearlet generator y/™" by 

V>r(^)-^/-;/2(<^)r(<^), 

where Pe{^) — f (2^+'^i,^2) for ^ > and the trigonometric polynomial P is a 2D 
fan filter (c.f. 1 12 1). For an illustration of P we refer to Fig.[TO|a). This in turn defines 
shearlets (/""m generated by non-separable generator functions for each scale 
index j > by setting 

where Mcj is a sampling matrix given by Mcj = diag(cj,C2) and c\ and are sam- 
pling constants for translation. 

One major advantage of these shearlets WjTm '■^^ ^^^^ ^^^^ ^ filter en- 
ables refinement of the directional selectivity in frequency domain at each scale. 



Fig. 10 a)-(b) show the refined essential support of y/]Tm compared to shearlets 



Yj,k,m arising from a separable generator as in Subsection 3.2.1 



Digital Shearlet Transforms 



39 



(a) 





Fig. 10 (a) Magnitude response of 2D fan filter. (b)Non-separable shearlet V'/°" 
shearlet V;\/t.m- 



(c)Separable 



3.3.2 Algorithmic Realization 



Next, our aim is to derive a digital formulation of the shearlet coefficients (/, i/J^ 
for a function / as given in (l39|. We will only discuss the case of shearlet coefficients 



,non 



associated with and 5^; the same procedure can be appUed for and except 
for switching the order of variables xi and X2- 



In Subsection 3.2.3 we discretized a sheared function f{S2-j/2j^-) using the dig- 
ital shear operator 'Sj-j/zj- defined in (|46|. In this implementation, we walk a 
different path. We digitize the shearlets W'jTm(') — V]Tmi^2-j/^-k') combining 
multiresolution analysis and digital shear operator S'^_j^2^. ^'^ digitize the wavelet 
o"m ^^'^ shear operator 52-j/2^, respectively. This yields digitized shearlet fil- 
ters of the form 



where Wj is the 2D separable wavelet filter defined by Wj{ni,n2) = gj-j{ni)- 
hj-j/iini) and pj-niin) are the Fourier coefficients of the 2D fan filter Pj-ji2- 
The DNST associated with the non-separable shearlet generators 1/"°" is then given 
by 

DNSTj^,{fj){n) = (/y* v7d,)(2^-.'c{«i,2^-^/24«2), for/y G ^^(Z^). 

We remark that the discrete shearlet filters Xj/j^. are computed by using a similar 
ideas as in Subsection |3.2.1| As before, those filter coefficients can be precomputed 
to avoid additional computational effort. 



Further notice that by setting cj = 2' 



and 



2'/2--', the DNST simply be- 



comes a 2D convolution. Thus, in this case, DNST is shear invariant. 



3.3.3 Inverse DNST 



In case that Cj = 2-' and = 2'/^ ■', the dual shearlet filters xj/j^. can be easily 
computed by deconvolution, and we obtain the reconstruction formula 
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Thus, no iterative methods are required for the inverse DNST. 

The frequency response of a discrete shearlet filter Xj/jj^ and its dual xj/j/^ is illus- 



trated in Fig. 1 1 We observe that primal and dual shearlet filters behave similarly in 
the sense that both of filters are very well localized in frequency. 



Fig. 11 Magnitude response of shearlet filter and its dual filter xjfj^- 



4 Framework for Quantifying Performance 

We next present the framework for quantifying performance of implementations of 
directional transforms, which was originally introduced in ll20l[T4l . This set of test 
measures was designed to analyze particular well-understood properties of a given 
algorithm, which in this case are the desiderata proposed at the beginning of this 
chapter. This framework was moreover introduced to serve as a tool for tuning the 
parameters of an algorithm in a rational way and as an objective common ground 
for comparison of different algorithms. The performance of the three algorithms 
FDST, DSST, and DNST will then be tested with respect to those measures. This 
will give us insight into their behavior with respect to the analyzed characteristics, 
and also allow a comparison. However, the test values of these three algorithms will 
also show the delicateness of designing such a testing framework in a fair manner, 
since due to the complexity of algorithmic realizations it is highly difficile to do 
each aspect of an algorithm justice. It should though be emphasized that - apart 
from being able to rationally tune parameters - such a framework of quantifying 
performance is essential for an objective judgement of algorithms. The codes of all 
measures are available in lShearLab 

In the following, S shall denote the transform under consideration, S* its adjoint, 
and, if iterative reconstruction is tested, G^J shall stand for the solution of the matrix 
problem AI ~ J using the conjugate gradient method with residual error set to be 
10^^. Some measures apply specifically to transforms utilizing the pseudo-polar 
grid, for which purpose we introduce the notation P for the pseudo-polar Fourier 
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transform, w shall denote the weighting applied to the values on the pseudo-polar 
grid, and W shall be the windowing with additional 2D iFFT. 



4.1 Algebraic Exactness 

We require the transform to be the precise implementation of a theory for digital data 
on a pseudo-polar grid. In addition, to ensure numerical accuracy, we provide the 
following measure. This measure is designed for transforms utilizing the pseudo- 
polar grid. 

Measure 1 Generate a sequence of 5 (of course, one can choose any rea- 
sonable integer other than 5) random images /i,...,/5 on a pseudo-polar 
grid for N = 512 and R — 8 with standard normally distributed entries. Our 
quality measure will then be the Monte Carlo estimate for the operator norm 
\\W*W-Id\\op given by 

\\W*WI,-Ii\\2 

Maig = max . 

'=l,-,5 ||/,i|2 



This measure applies to the FDST - not to the DSST or DNST - , for which we 
obtain 

Maif, = 6.6E - 16. 

This confirms that the windowing in the FDST is indeed up to machine precision a 
Parseval frame, which was already theoretically confirmed by Theorem |2] 

4.2 Isometry of Pseudo-Polar Transform 

We next test the pseudo-polar transform itself which might be used in the algorithm 
under consideration. For this, we will provide three different measures, each being 
designed to test a different aspect. 

Measure 2 

• Closeness to isometry. Generate a sequence of 5 random images Ii,. . . ,1$ 
of size 512x512 with standard uniformly distributed entries. Our quality 
measure will then be the Monte Carlo estimate for the operator norm 
\\P*wP — Id\\ap given by 
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\\P*wPIi^U\\2 

Misatrn max T—- . 

'=l,-,5 ||/,-||2 

• Quality of preconditioning. Our quality measure will be the spread of the 
eigenvalues of the Gram operator P*wP given by 

_ ^nax(f*wf ) 

• Invertibility. Our quality measure will be the Monte Carlo estimate for the 
invertibility of the operator \/wP using conjugate gradient method G 
{residual error is set to be 10^^, here GaJ means solving matrix problem 
AI = J using conjugate gradient method) given by 



\G^py^PIi-Ii\\l 



Misomi = max 



This measure applies to the FDST - not to the DSST or DNST - , for which we 
obtain the following numerical results, see Table |5] 



Table 5 The numerical results for the test on isometry of the pseudo-polar transform. 











FDST 


9.3E-4 


1.834 


3.3E-7 



The slight isometry deficiency of Misomi ~ 9.9E-4 mainly results from the isom- 
etry deficiency of the weighting. However, for practical purposes this transform can 
be still considered to be an isometry allowing the utihzation of the adjoint as inverse 
transform. 



4.3 Parseval Frame Property 

We now test the overall frame behavior of the system defined by the transform. 
These measures now apply to more than pseudo-polar based transforms, in particu- 
lar, to FDST, DSST, and DNST. 

Measure 3 Generate a sequence of 5 random images Ii, . . . ,1^ of size 512 x 
512 with standard uniformly distributed entries. Our quality measure will then 
two-fold: 
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• Adjoint transform. The measure will be the Monte Carlo estimate for the 
operator norm \\S* S — Id\\op given by 

\\S*SIi-Ii\\2 

Mnghti = max . 

'=1 5 ||y/||2 

• Iterative reconstruction. Using conjugate gradient G ^,p, our measure 
will be given by 

"""^"'^ = ,.^"5 Ml ■ 



The following table. Table |6] presents the performance of FDST and DNST with 
respect to these quantitative measures. 



Table 6 The numerical results for the test on Parseval property. 







Mtight2 


FDST 


9.9E-4 


3.8E-7 


DSST 


1.9920 


1.2E-7 


DNST 


0.1829 


5.8E-16 (with dual filters) 



The transform FDST is nearly tight as indicated by the measures Mtighti = 9.9E- 
4, i.e., the chosen weights force the PPFT to be sufficiently close to an isom- 
etry for most practical purposes. If a higher accurate reconstruction is required, 
Mtight^ — 3.8E-7 indicates that this can be achieved by the CG method. As expected, 
Mtight, = 1.9920 shows that DSST is not tight. Nevertheless, the CG method pro- 
vides with Mtightj = 1 ■2E-7 a highly accurate approximation of its inverse. DNST is 
much closer to being tight than DSST (see Mtight, = 0.1829). We remark that this 
transform - as discussed - does not require the CG method for reconstruction. The 
value 5.8-16 was derived by using the dual shearlet filters, which show superior 
behavior. 



4.4 Space-Frequency-Localization 



The next measure is designed to test the degree to which the analyzing elements, 
here phrased in terms of shearlets but can be extended to other analyzing elements, 
are space-frequency localized. 
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Measure 4 Let I be a shearlet in a 512 x 512 image centered at the origin 
(257,257) with slope of scale 4, i.e., O4QQ + C74QQ. Our quality measure 
will be four-fold: 

• Decay in spatial domain. We compute the decay rates di,... ,^512 along 
lines parallel to the y-axis starting from the line [257, : ] and the de- 
cay rates dsn, ■ ■ •. ^^1024 with x and y interchanged. By decay rate, for 
instance, for the line [257 : 512, 1], we first compute the smallest mono- 
tone majorant M{x, 1), x = 257, ... ,512 - note that we could also choose 
an average amplitude here or a different 'envelope' - for the curve 
|/(x, 1)|, X = 257, . . . ,512. Then the decay rate is defined to be the aver- 
age slope of the line, which is a least square fit to the curve log(M(x, 1)), 
X = 257, ... ,512. Based on these decay rates, we choose our measure to 
be the average of the decay rates 

Mdecayi = ^ d-i- 

lUZt ,=1 1024 



Decay in frequency domain. Here we intend to check whether the Fourier 
transform of I is compactly supported and also the decay. For this, let I 
be the 2D-FFT of I and compute the decay rates d,, i = 1, . . . , 1024 as 
before. Then we define the following two measures: 
- Compactiy supportedness. 



- Decay rate. 



_ max|»|,|v|<3 \I{u,v) 
max„,v |^(«,v)| 



Smoothness in spatial domain. We will measure smoothness by the aver- 
age of local Holder regularity. For each (mq, vq), we compute M{u,v) = 
|/(m, v) — /(mo, vo)|, < max{|M — Mo|, |v — vo|} < 4. Then the local Holder 
regularity O£„p^vo the least square fit to the curve log(|M(M, v)|). Then 
our smoothness measure is given by 

1 ^ 

^smoothi — cTo2 £^^tV 
•'^^ u,v 

Smoothness in frequency domain. We compute the smoothness now for I, 
the 2D-FFTofI to obtain the new a^^v ctnd define our measure to be 

1 ^ 

^smooth2 — cTo2 w ' 
•'^^ «,v 
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Let us now analyze the space-frequency localization of the shearlets utilized in 
FDST, DSST and DNST by these measures. The numerical results are presented in 
Table El 



Table 7 The numerical results for the test on space-frequency localization. 







Msupp 


^decay2 


^smooth I 


^smooth 2 


FDST 


-1.920 


5.5E-5 


-3.257 


1.319 


0.734 


DSST 


— oo 


8.6E-3 


-1.195 


0.012 


0.954 


DNST 


— oo 


2.0E-3 


-0.716 


0.188 


0.949 



The shearlet elements associated with FDST are band-limited and those asso- 
ciated with DSST and DNST are compactly supported, which is clearly indicated 
by the values derived for Mciecayi^ Ms„pp, and Mdecay2- It would be expected that 
Mdecay2 — ~°° ^r FDST due to the band-limitedness of the associated shearlets. 
The shearlet elements are however defined by their Fourier transform on a pseudo- 
polar grid, whereas the measure Mii^cuyi is taken after applying the 2D-FFT to the 
shearlets resulting in data on a cartesian grid, in particular, yielding a non-precisely 
compactly supported function. 

The test values for Mv,„o„,/,j and Msmooth2 show that the associated shearlets are 
more smooth in spatial domain for FDST than for DSST and DNST, with the re- 
versed situation in frequency domain. 



4.5 True Shear Invariance 

Shearing naturally occurs in digital imaging, and it can - in contrast to rotation - be 
precisely realized in the digital domain. Moreover, for the shearlet transform, shear 
invariance can be proven and the theory implies 

We therefore expect to see this or a to the specific directional transform adapted be- 
havior. The degree to which this goal is reached is tested by the following measure. 

Measure 5 Let I be an 256 x 256 image with an edge through the origin 
(129, 129) of slope 0. Given —l<s<\, generates an image Is '■= I{Ss-) and 
let Sj be the set of all possible scales j such that l-^s € Z. Our quality measure 
will then be the curve 
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M shear J ^ max r-r , scalejGSj, 

-2i<k,k+2Js<2i II-/ II 2 

where Cj± is the shearlet coefficients at scale j and shear k. 

We present our results in Table |8] 
Table 8 The numerical results for the test on shear invariance. 







^ shear. 2 




^shear.A 


FDST 


1.6E-5 


1.8E-4 


0.002 


0.003 



This table shows that the FDST is indeed almost shear invariant. A closer inspec- 
tion shows that Mshear,\ ™d Mshear.2 ^re relatively small compared to the measure- 
ments with respect to finer scales Mshear.i ™d M shear A- The reason for this is the 
aliasing effect which shifts some energy to the high frequency part near the bound- 
ary away from the edge in the frequency domain. 

We did not test DSST and DNST with respect to this measure, since these trans- 
forms show a different - not included in this Measure 5 - type of shear invariance 
behavior. 



4.6 Speed 



Speed is one of the most fundamental properties of each algorithm to analyze. Here, 
we test the speed up to a size of = 512 which regard as sufficient to computing 
the complexity. 



Measure 6 Generate a sequence of 5 random images If, / = 5, . . . ,9 of size 
2' X 2' with standard normally distributed entries. Let Sj be the speed of the 
shearlet transform S applied to !(. Our hypothesis is that the speed behaves 
like Si = c ■ {i^'Y; 2^' being the size of the input. Let now da be the average 
slope of the line, which is a least square fit to the curve i i-^ log(i, ). Let also 
fi be the 2fft applied to It, i = 5,..., 9. Our quality measure will then be three- 
fold: 

• Complexity. 

^^P^^^' "21og2- 

• The Constant 
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If. S, 

Comparison with 2D-FFT. 

,=5 



Table |9] presents the resuhs of testing FDST, DSST and DNST with respect to 
these speed measures. 



Table 9 The numerical results for the test on speed. 





^ Speed 1 




^speedi, 


FDST 


1.156 


9.3E-6 


280.560 


DSST 


0.821 


4.5E-3 


88.700 


DNST 


1.081 


9.9E-8 


40.519 



To interpret these results correctly, we remark that the DNST was tested only with 
test images /, for / = 7, ... ,9, since it can not be implemented for small size images. 
Interestingly, the results also show that the 2D FFT based convolution makes DNST 
comparable to DSST with respect to these speed measures, although it is much 
more redundant than DSST. Finally, the results show that FDST is comparable with 
both DSST and DNST with respect to complexity measure Mspeedp From this, it 
is conceivable to assume that FDST is highly comparable with respect to speed 
for large scale computations. The larger value Mspeed^ = 280.560 appears due to 
the fact that the FDST employs fractional Fourier transforms on an oversampled 
pseudo-polar grid of size. 



4.7 Geometric Exactness 



One major advantage of directional transforms is their sensitivity with respect to 
geometric features alongside with their ability to sparsely approximate those (cf. 
Chapter Q). This measure is designed to analyze this property. 

Measure 7 Let Ii,. . . ,1^ be 256 x 256 images of an edge through the origin 
(129,129) and of slope [—1,-0.5,0,0.5,1] and the transpose of the middle 
three, and let cij be the associated shearlet coefficients for image li and scale 
j. Our quality measure will two-fold: 
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Decay of significant coefficients. Consider the curve 
1 ^ 

- ^ max \ci j(of analyzing elements aligned with the line)\, scale 7, 

let d be the average slope of the line, which is a least square fit to log of 
this curve, and define 

Mgeoi = d. 



• Decay of insignificant coefficients. Consider the curve 
1 ^ 

-^max\cjj(of all other analyzing elements)], scale j, 
° 1=1 

let d be the average slope of the line, which is a least square fit to log of 
this curve, and define 

Mgeo2 = d. 



Table[TO]shows the numerical test results for FDST, DSST, and DNST. 
Table 10 The numerical results for the test on geometric exactness. 









FDST 


-1.358 


-2.032 


DSST 


-0.002 


-0.030 


DNST 


-0.019 


-0.342 



As expected, the decay rate of the insignificant shearlet coefficients of FDST, 
i.e., the ones not aligned with the line singularity, measured by M^go^ ~ -2.032 is 
much larger than the decay rate of the significant shearlet coefficients measured by 
Mgecu ~ -1.358. Notice that this difference is even more significant in the case of 
the DSST and DNST. 



4.8 Robustness 

To analyze robustness of an algorithm, we choose thresholding as the most common 
impact on a sequence of transform coefficients. 



Measure 8 Let I be the regular sampling of a Gaussian function with mean 
and variance 256 on {—128, 127}^ generating an 256 x 256-image. 
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• Thresholding 1. Our first quality measure will be the curve 

_ ||G^^W*thres^p,5/-/||2 

'^'hres^^Pl - |j/|j2 

where thresi pj discards 100 • (1 — 2^''' ) percent of the coefficients (p\ = 
[2:2: 10]j. ' 

• Thresholding 2. Our second quality measure will be the curve 

_ HGy^fW*thres2,p,5/~/||2 

Mr/,rei2,P2 ~ ||/||2 

where thres2,p2 ^^t^ '^(/Z those coefficients to zero with absolute values be- 
low the threshold m(l — 2^''^) with m being the maximal absolute value 
of all coefficients. (p2 = [0.001 : 0.01 : 0.041] j 



Table 



11 



shows that even if we discard 100(1 - 2"'°) - 99.9% of the FDST co- 
efficients, the original image is still well approximated by the reconstructed image. 
Thus the number of the significant coefficients is relatively small compared to the 
total number of shearlet coefficients. From Table [12] we note that knowledge of 
the shearlet coefficients with absolute value greater than m(l — 1 /2'' '"^^)(^ 0.1% of 
coefficients) is sufficient for precise reconstruction. 

DNST shows a similar behavior with worse values for relatively large pi. It 
should be however emphasized that firstly, the redundancy of DNST used in this 
test is 25 and this is lower than the redundancy of FDST, which is about 71. This 
effect can be more strongly seen by the test results of DSST whose redundancy with 
4 even much smaller. Secondly, a significant part of the low frequency coefficients 
in both DSST and DNST will be removed by a relatively large threshold, since the 
ratio between the number of the low frequency coefficients and the total number of 
coefficients is much higher than FDST. This prohibits a similarly good reconstruc- 
tion of a Gaussian function. 

This test in particular shows the delicateness of comparing different algorithms 
by merely looking at the test values without a rational interpretation; in this case, 
without considering the redundancy and the ratio between the number of the low 
frequency coefficients and the total number of coefficients. 



Table 11 The numerical results for M,itresj 



P\ 


2 


4 


6 


8 


10 


FDST 


1.5E-08 


7.2E-08 


2.5E-05 


0.001 


0.007 


DSST 


0.02961 


0.02961 


0.02961 


0.0296 


0.0331 


DNST 


5.2E-10 


1.2E-04 


0.00391 


0.0124 


0.0396 
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Table 12 The numerical results for M,f,res2 



P2 0.001 0.011 0.021 0.031 0.041 



FDST 0.005 0.039 0.078 0.113 0.154 

DSST 0.030 0.036 0.046 0.056 0.072 

DNST 0.002 0.018 0.035 0.055 0.076 
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